xxxxxxxxxxbegin using Plots using PlutoUI using StatisticsendDeriving a SIR epidemiological model
In this notebook we will go from a discrete epidemiological model to its continuous counterpart. We will go from a (discrete) system of difference equiations to a system of (continuous) differential equations. But first: SIR.
Susceptible, Infected, Recovered (SIR)
In this workhorse of epidemiology we model agents as belonging to a certain class, or compartment. Either they are
susceptible to get infected,
infected, or
removed/recovered
The models you may have heard about during the COVID19 pandemic are more sophisticated versions of what we will see here, but they share some important common features. At the end of this notebook, we will want to generate the following plot. It shows the number of agents belonging to each class at a certain point in time after the outbreak of the disease:
xxxxxxxxxxdiscrete_time_SIR_plotWe will make this very simple. Instead of trying to model in great detail which type of agent (old, young, where, what time of the day, doing which activity etc) is likely to get infected, we will try to condense the overall probability of infection into a single number. Let's start to model the recovery from an infection
we have
people infected at time . If each has probability to recover each day, how many are still infected at day number ?
For any given person in
bernoulli (generic function with 1 method)xxxxxxxxxxbernoulli(p) = rand() < p 200xxxxxxxxxxN = 200step! (generic function with 1 method)xxxxxxxxxxfunction step!(infectious, p) for i in 1:length(infectious) if infectious[i] && bernoulli(p) infectious[i] = false end end return infectiousend100xxxxxxxxxxT = 100true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
false
true
true
true
true
false
true
true
true
false
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
false
false
false
true
true
false
true
true
true
false
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
false
false
false
true
true
false
false
true
true
false
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
false
false
false
true
true
false
false
true
true
false
true
true
true
true
true
true
false
true
true
true
true
true
true
true
true
true
true
true
true
true
false
false
false
true
true
false
false
true
true
false
true
true
true
true
true
true
false
true
true
false
true
true
true
true
true
true
true
true
true
true
false
false
false
true
true
false
false
true
true
false
true
true
true
true
true
true
false
true
true
false
true
true
false
true
true
true
true
true
true
true
false
false
false
true
true
false
false
true
true
false
true
true
true
true
true
true
false
true
true
false
true
false
false
true
true
true
false
true
true
true
false
false
false
true
true
false
false
false
true
false
true
true
true
true
true
true
false
true
true
false
true
false
false
true
true
true
false
true
true
true
false
false
false
true
true
false
false
false
true
false
true
true
true
true
true
true
false
true
true
false
true
false
false
false
true
true
false
true
true
true
false
false
false
true
true
false
false
false
true
false
true
false
true
true
true
true
false
false
true
false
true
false
false
false
true
true
false
true
true
true
false
false
false
true
true
false
false
false
true
false
true
false
true
true
true
true
false
false
true
false
true
false
false
false
false
true
false
true
false
true
false
false
false
true
true
false
false
false
true
false
true
false
true
true
true
true
false
false
true
false
true
false
false
false
false
true
false
true
false
true
false
false
false
true
true
false
false
false
true
false
true
false
true
true
true
true
false
false
true
false
true
false
false
false
false
true
false
true
false
true
false
false
false
true
true
false
false
false
true
false
true
false
true
true
true
true
false
false
true
false
true
false
false
false
false
true
false
true
false
true
false
false
false
true
true
false
false
false
true
false
true
false
true
true
true
true
false
false
true
false
true
false
false
false
false
true
false
true
false
true
false
false
false
true
true
false
false
false
true
false
true
false
true
true
true
true
false
false
true
false
true
false
false
false
false
true
false
true
false
true
false
false
false
true
true
false
false
false
false
false
true
false
true
true
true
true
false
false
true
false
true
false
false
false
false
true
false
true
false
true
false
false
false
true
true
false
false
false
false
false
true
false
false
true
true
true
false
false
true
false
true
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
true
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
true
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
true
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
true
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
xxxxxxxxxx# try this out: all are infected on day one.# let's record how many are in each state I or R after i daysbegin infected = trues(N) # all infected results = [copy(step!(infected, 0.05)) for i in 1:T] # we `copy` here because we keep modifying the same array throughout pushfirst!(results, trues(N))endxxxxxxxxxx i Slider(1:T, show_value=true)ok great, now we know how many are count how many infected there are each day!
simulate_recovery (generic function with 1 method)xxxxxxxxxxfunction simulate_recovery(p, T) infectious = trues(N) num_infectious = [N] for t in 1:T step!(infectious, p) push!(num_infectious, count(infectious)) end return num_infectiousendxxxxxxxxxxbegin pp = 0.05 plot(simulate_recovery(pp, T), label="run 1", alpha=0.5, lw=2, m=:o) plot!(simulate_recovery(pp, T), label="run 2", alpha=0.5, lw=2, m=:o) xlabel!("time t") ylabel!("number infectious")endxxxxxxxxxxall_data = [simulate_recovery(pp, T) for i in 1:30]; # do 30 runsxxxxxxxxxxbegin plot(replace.(all_data, 0.0 => NaN), yscale=:log10, alpha=0.3, leg=false, m=:o, ms=1, size=(500, 400)) plot!(mean(all_data), yscale=:log10, lw=3, c=:red, m=:o, label="mean", alpha=0.5) xlabel!("time t") ylabel!("number infectious")endDeterministic dynamics for the mean: Intuitive derivation
The mean seems to behave in a rather predictable way over time. Can we derive this?
Let
So
or
At time
We can rearrange and solve the recurrence relation:
and hence solve the recurrence relation:
Let's compare the exact and numerical results:
200.0
190.0
180.5
171.475
162.901
154.756
147.018
139.667
132.684
126.05
119.747
113.76
108.072
102.668
97.535
92.6582
88.0253
83.6241
79.4429
75.4707
1.87879
1.78485
1.69561
1.61083
1.53029
1.45377
1.38108
1.31203
1.24643
1.18411
xxxxxxxxxxbegin plot(mean(all_data), m=:o, label="numerical mean") plot!(exact, lw=3, label="analytical mean")end They agree well, as they should. The agreement is expected to be better (i.e. the fluctuations smaller) for a larger population.
Continuous time
If we look at the graph of the mean as a function of time, it seems to follow a smooth curve. Indeed it makes sense to ask not only how many people have recovered each day, but to aim for finer granularity.
Suppose we instead increment time in steps of
Then we will need to adjust the probability of recovery in each time step. It turns out that to make sense in the limit
where
We get
Dividing by
We recognise the left-hand side as the definition of the derivative when
That is, we obtain an ordinary differential equation that gives the solution implicitly. We know this has a general solution of the form
and together with an initial value of
Starting at time zero (
SIR model
Now let's extend the procedure to the full SIR model,
xxxxxxxxxxmd"""Now let's extend the procedure to the full SIR model, $S \to I \to R$. Since we already know how to deal with recovery, consider just the SI model, where susceptible agents are infected via contact, with a certain probability."""Let's denote by
On average, in each sweep each infectious individual has the chance to interact with one other individual. That individual is chosen uniformly at random from the total population of size
Hence the change in the number of infectious people after that step is:
It is useful to normalize by
Including recovery with probability
xxxxxxxxxxmd"""$$\begin{align}s_{t+1} &= s_t - b \, s_t \, i_t \\i_{t+1} &= i_t + b \, s_t \, i_t - c \, i_t\\r_{t+1} &= r_t + c \, i_t\end{align}$$"""And again we can allow the processes to occur in steps of length
xxxxxxxxxxmd"""$$\begin{align}\textstyle \frac{ds(t)}{dt} &= -\beta \, s(t) \, i(t) \\\textstyle \frac{di(t)}{dt} &= +\beta \, s(t) \, i(t) &- \gamma \, i(t)\\\textstyle \frac{dr(t)}{dt} &= &+ \gamma \, i(t)\end{align}$$"""Note that no analytical solutions of these (simple) nonlinear ODEs are known as a function of time!
Below is a simulation of the discrete-time model. Note that the simplest numerical method to solve (approximately) the system of ODEs, the Euler method, basically reduces to solving the discrete-time model! A whole suite of more advanced ODE solvers is provided in the Julia DiffEq ecosystem. We will in another notebook introduce that package, together with a continuous version of the SIR model.
xxxxxxxxxxbegin ts = 1:length(SIR) discrete_time_SIR_plot = plot(ts, [x.s for x in SIR], m=:o, label="S", alpha=0.2, linecolor=:blue, leg=:right, size=(400, 300)) plot!(ts, [x.i for x in SIR], m=:o, label="I", alpha=0.2) plot!(ts, [x.r for x in SIR], m=:o, label="R", alpha=0.2) xlims!(0, 500)end0xxxxxxxxxxbegin NN = 100 SS = NN - 1 II = 1 RR = 0end0.99
0.01
0.0
xxxxxxxxxxss, ii, rr = SS/NN, II/NN, RR/NN0.1
0.01
xxxxxxxxxxp_infection, p_recovery = 0.1, 0.01discrete_SIR (generic function with 2 methods)xxxxxxxxxxfunction discrete_SIR(s0, i0, r0, T=1000) s, i, r = s0, i0, r0 results = [(s=s, i=i, r=r)] for t in 1:T Δi = p_infection * s * i Δr = p_recovery * i s_new = s - Δi i_new = i + Δi - Δr r_new = r + Δr push!(results, (s=s_new, i=i_new, r=r_new)) s, i, r = s_new, i_new, r_new end return resultsend0.99
0.01
0.0
0.98901
0.01089
0.0001
0.987933
0.0118581
0.0002089
0.986761
0.0129111
0.000327481
0.985487
0.014056
0.000456592
0.984102
0.0153006
0.000597151
0.982597
0.0166533
0.000750157
0.98096
0.0181231
0.000916691
0.979182
0.0197197
0.00109792
0.977251
0.0214534
0.00129512
0.975155
0.0233354
0.00150965
0.972879
0.0253777
0.00174301
0.97041
0.0275928
0.00199678
0.967733
0.0299945
0.00227271
0.96483
0.0325973
0.00257266
0.961685
0.0354164
0.00289863
0.958279
0.0384681
0.00325279
0.954593
0.0417698
0.00363748
0.950605
0.0453394
0.00405517
0.946295
0.049196
0.00450857
3.65235e-5
8.28801e-5
0.999881
3.65232e-5
8.20516e-5
0.999881
3.65229e-5
8.12314e-5
0.999882
3.65226e-5
8.04194e-5
0.999883
3.65223e-5
7.96155e-5
0.999884
3.65221e-5
7.88196e-5
0.999885
3.65218e-5
7.80317e-5
0.999885
3.65215e-5
7.72517e-5
0.999886
3.65212e-5
7.64794e-5
0.999887
3.65209e-5
7.57149e-5
0.999888
xxxxxxxxxxSIR = discrete_SIR(ss, ii, rr)